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Abstract 

In this paper we present a new multi-asset pricing model, which is built upon newly developed 
families of solvable multi-parameter single-asset diffusions with a nonlinear smile-shaped volatil- 
ity and an afHne drift. Our multi-asset pricing model arises by employing copula methods. In 
particular, all discounted single-asset price processes are modeled as martingale diffusions under a 
risk-neutral measure. The price processes are so-called UOU diffusions and they are each generated 
by combining a variable (Ito) transformation with a measure change performed on an underlying 
Ornstein-Uhlenbeck (Gaussian) process. Consequently, we exploit the use of a normal bridge copula 
for coupling the single-asset dynamics while reducing the distribution of the multi-asset price pro- 
cess to a multivariate normal distribution. Such an approach allows us to simulate multidimensional 
price paths in a precise and fast manner and hence to price path-dependent financial derivatives 
such as Asian-style and Bermudan options using the Monte Carlo method. We also demonstrate 
how to successfully calibrate our multi-asset pricing model by fitting respective equity option and 
asset market prices to the single-asset models and their return correlations (i.e. the copula function) 
using the least-square and maximum-likelihood estimation methods. 

Introduction 

Many quantitative finance applications require a multi-asset pricing model vi^ith dependencies between 
the single-asset price components. Compared to the variety of univariate asset price models, the pool of 
multi-asset pricing models is not so extensive. Most of multivariate models are based on multidimensional 
geometric Brownian motion with the possible inclusion of a jump process. In this paper we develop and 
explore a new multi-asset arbitrage-free pricing model based on a special family of nonlinear diffusions. 
The development of efficient computational methods for pricing multi-asset equity derivatives under 
such a model and the calibration of the multi-asset model to both standard equity option data as well 
as historical equity prices are the objectives of the current paper. 

Here, we specialize on option pricing applications under so-called UOU diffusion models, which are 
obtained by transforming an underlying Ornstein-Uhlenbeck diffusion process via the use of a diffusion 
canonical transformation method (see [21 [31 |SJ [5] and references therein). For all choices of model 
parameters, all discounted (single-asset) price processes UOU are conservative martingales under a risk- 
neutral measure. Since the univariate diffusions are solvable, the single-asset risk-neutral transition 
probability density function is given in analytically closed form. Moreover, implied volatility surfaces for 
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this highly nonlinear asset price model exhibit a wide range of pronounced smiles and skews of the type 
observed in the option markets. The main relevant features of the univariate UOU model are summarized 
in Section [H 

To construct a multivariate probability distribution, one can use a copula function that allows us 
to couple univariate distribution functions. Sampling from the obtained joint multivariate distribution 
function thereby reduces to sampling from the copula function and from the univariate distributions. 
Therefore, the copula method allows us to construct the joint distribution and density functions as well 
as to obtain an exact path sampling algorithm. 

The main computational disadvantage of such an approach is the calculation of inverses of the distri- 
bution functions. This operation can be a rather time-consuming computational problem for a compli- 
cated multi-parameter distribution. Nevertheless, such a drawback can be significantly improved if the 



copula function and univariate distributions have a similar structure. As is shown in Subsection 1.3 the 



bridge probability density function (conditional on the values of the process at the endpoints of a time 
interval) of a UOU diffusion is reduced to a normal density. Hence it is natural to couple univariate UOU 
bridges using a Gaussian copula. Based on this idea, in Section [2j we construct a two-step algorithm for 
the exact path simulation of the multidimensional (nonlinear) UOU process in the risk- neutral measure. 
Firstly, we apply a usual copula method for sampling the multi-asset process at the terminal (maturity) 
time. Secondly, we use a bridge sampling along with a multivariate normal distribution to model the 
process at any intermediate time. 

In Section [3] we demonstrate the calibration of the univariate and multivariate models to historical 
asset and equity option prices. The calibration process has two stages. First, we calibrate all univariate 
(marginal) asset price models independently of each other. Using the least-square method the models 
can be fitted to standard European option prices. Alternatively, the maximum likelihood estimation 
(MLE) allows the models to be fitted to historical asset prices. Second, we fit the copula function to 
historical observations. Since, our model assumes a normal copula, we need to find a best-fitted normal 
correlation matrix. 

In Section [4j we give computational applications of the model to pricing multi-asset path-dependent 
Asian-style and Bermudan options. In pricing Bermudan options we use a regression-based Monte Carlo 
method. 

In summary, the main results of our paper include: the construction a new family of multivariate 
models for which marginal processes are local volatility smile diffusions; the development of calibration 
schemes for the single-asset and multi- asset pricing UOU diffusion models based on the least-square and 
MLE methods; the construction of an exact multivariate path simulation method that can be used for 
Monte Carlo pricing of generally path-dependent European-style and American-style options. 



1 Ornstein-Uhlenbeck Family of Univariate State-Dependent 
Volatility Diffusion Models 

1.1 Diffusion Canonical Transformation 

The diffusion canonical transformation method, first presented in \T and then further developed in 
EJ H], leads to various families of solvable one-dimensional diffusions with a nonlinear diffusion 
coefficient function and an affine drift. In this paper, we consider the UOU family, which is based on 
a regular Ornstein-Uhlenbeck process {Xt)t>Q E I = R. The regular Ornstein-Uhlenbeck process is 
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defined by the infinitesimal generator 

£f{x) = ^iy^r{x)-Xxf'{x), xeR, (1.1) 

where A > and v > are constants. Both left and right boundaries / = — oo and r = oo of the 
state space I are non-attracting natural for all choices of parameters. The transition probability density 
hmction (PDF) is 

, , I K / k{x — xoe~^*)'^\ ^, 

^^(''"°'") = VMl-e-A.)^^P(,- 2(1 -e-A*) ) ' ^'-'^ 

where we define k = ^ > 0. 

Let p > be a strictly positive constant. Then, two fundamental solutions oi Cip{x) — sip{x), x €l, 
are given by 

ip-{x) =expi—j D^^{x^/^) and tp+ (x) = tp' (-x), (1.3) 

where v = p/X > and D_^{-) is Whittaker's parabolic cylinder function (see [T]). The solutions 
t/jp (x) and iPp{x) are linearly independent and respectively increasing and decreasing positive functions 
of a; G E. 

We now construct another diffusion process {St)t>a G 2? = = (0,oo) by applying a diffusion 
canonical transformation to the underlying Ornstein-Uhlcnbeck diffusion. This process obeys the stochas- 
tic differential equation (SDE) 

dSt = rStdt + aiSt)dWu St=o = So , (1.4) 

where r is a constant and <y{S) is a nonlinear diffusion coefficient function. Here {Wt)t>Q is a standard 
Brownian motion in some appropriate probability measure P. 

The initial step of the transformation is to apply a Doob's /i-transform or p-excessive transform (see 
[4]) to (Xt). The resulting diffusion process (Xj''-')^^^^ is defined by the following infinitesimal generator: 

C^'^fix) = I'^'fix) + Xp{x)f{x), xel, (1.5) 

where Xp{x) = —Xx + i''^Up{x)/up{x) for any p > 0. The strictly positive generating function is given 

by Up{x) = qiip'^{x) + q2ipp{x) with parameters qi,q2 > 0, qi + q2 > 0. The process (X^^^) has the 
transition PDF 

py{t;xo,x) = e~''*^r4 — vPx {t]Xo,x) , x,xo el, t > 0. (1.6) 

^PV^O ) 

The final step of the transformation is a change of variabl e (se e [S] for a more general discussion). 
We define a new process St = f^X^^'^), t>0, that solves SDE (1.4) by finding a strictly monotonic map 
that solves C^p^F{x) = rF{x), for constant r. Then a{F{x)) ~ i'\F'{x)\ or equivalently a{s) ~ i^/|X'(s)|, 
where X = defines the unique inverse map. The transition PDF "ps for the process (St) follows from 
that for the underlying process [Xt): 

I \ (p) I \ ^ Up(x) _ f , , /-, -\ 

Ps{t]S(i,s) ^ —^py{t]Xo,x) ^ — re I" px{t;xQ,x) , (1.7) 
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X = X(s), xo = X(so). 

We now apply the above construction to a subfamily of diffusions with choice gi = 0,(72 = 1, i-e. 
with generating function Up(x) — Lp~{x). Letting r + p > and p > 0, we specifically consider 

Fa;) = c— =c — - — . (1.8) 

Ifp (x) D^^{Xy^K) 

This function maps a; e R onto s € (0, 00) and is monotonically increasing, where s = f{x) has the 
unique inverse relation x — X(s). This transformation leads to a family of diffusions {St) that is referred 
to as the unbounded Ornstein-Uhlcnbeck (UOU) family with the diffusion coefficient function given by 



where x ~ X(s). The volatility function a in (1.9 1 depends on several adjustable positive parameters 



such as c, I', A, p, and more generally r 7^ also enters as an additional parameter within the volatility 



specification. Notice that for the driftless case with r = formula (1.91 simplifies as follows: cr(F(x)) 
} , where uo > is a constant and s(x) — is the scale density for the X-diffusion. 
The following is an important statement for the purposes of risk-neutral pricing. 



Lemma 1.1 ([8 ). Consider a process {St)t>o G 1R,+ of the UOU family solving the SDE {1-4) with the 
ifficient a specified by (l.H^ and having transition PDF ps specified by equations { 1.1), ( 1 



usion coet 



1.9) and given by 



Ps{t;so,s)^ TTT^r^ jz—^, ^px{t;xo,x), t>0, s,so>0, (1.10) 

where W{x) = W[D^l,{x^/K), D_jj_r/\{—Xy/K)] > is the Wronskiaiyi and x — X(s), xq — X(so). 
Then (St) is a conservative process, i.e. P(5f G ]R+) = 1 for all t > O.Adoreover, the discounted asset 
price process {e~^^ St)t>o a true martingale. 

1.2 Pricing Vanilla Options 

According to Lemma |1.1| there exists an equivalent martingale measure under the UOU family with 



transition PDF ( 1.10 ) for any choice of the model parameters. Consider a standard European-style option 
defined by its payoff function A(S) at terminal price S = St and maturity (expiration) time T > to = 0. 
For example, a vanilla European call has the payoff function C^{S) = {S — K)^ = maxjS" — if, 0}, 
where if > is a strike price. The valuation of a standard European option is given by the conditional 
expectation under a risk-neutral probability measure P = Q: 

V{So,T) - e-'-^E«3[A(5T) | 5*^0 - ^0] = e-^^E*^ [A(F(A^^))) | X^"'^ = xq] . (1.11) 
This is reduced to the vahiation of a one-dimensional integral expressible as follows: 

V(So,T) = e-'-^ / ps{T;So,s)A{s)ds 
Jo 

Up{x) pxiT; xo, x)A{F{x)) dx , 



Up{xo) 



(1.12) 



-"-The Wronskain W of functions / and g is defined by W[f{x),g{x)] = f{x)g'{x) — f'(x)g{x). 
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Figure 1: Local volatility (or log-normal volatility) function (j{S)/S (left plot) and Black-Scholes impled 
volatility surface (right plot) for the UOU family. The UOU family is plotted using model parameters 
c = 100; p = 0.005, K = 5, u = 0.1 (thin line) and p = 0.02, n = 1, v = 0.5 (thick line). The implied 
volatility surface, with T as time to maturity and K as strike price, corresponds to the same choice 
of parameters as the local volatility plot drawn with the thin line. The interest rate r = 0.1 and spot 
5'o = 100 are used. 



where xq ~ X{So). 

1.3 The UOU Bridge 

Let us consider a bridge process on [ii,i2], < < t2, generated by a single asset price process St with 
St-i and St2 fixed at si and S2, respectively. The bridge density bs{t; s) = bs{ti,t2,t; si, S2, s) for St — s, 
ti < t < t2, conditional on St-^ = si and St2 — S2, is given by 

, f, X Ps{t~'ti;si,s)ps{t2~t;s,S2) 
bs{t;s) = — 



Psih -tl]Si,S2) 
V {t-h;xi, X{s))p^p {t2 - t; X(s) , X2) 



'^(^) P^P{t2-h;xuX2) 



(1.13) 



px{t - ti;xi,X{s))px{t2 - t;X{s),X2) _ v 



Pxih - h;xi,X2) 



a{s) 



bx{t;X{s)), 



where Xi = X(si), X2 = X(s2), and bxit;x) = bx{ti,t2,t;xi,X2,x) is the bridge PDF of the underlying 
X-diffusion conditional on the endpoint values Xt-^ — xi and Xt2 = X2- Since the underlying process 
{Xt) is a Gaussian process, the bridge process, generated by (S'f)t>o, with transition PDF (1.131 is just 



a nonlinear transformation of a Gaussian bridge, with respective path points mapped as Xt — X{St 
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Thus, the bridge PDF bs for (St) from the UOU family can be reduced to a Gaussian PDF for the 
underlying X-process. Indeed, from Eq. (1.2) we see that the bridge PDF bx of the Ornstein-Uhlenbeck 
diffusion is a normal density ^{x) 



1 e (a; fj.) /2b ^[^^^ mean u and variance b^ given by 



fe2 



xie 



AAi f^2\A2 



2AAi 



g2A(Ai+A2) 
^)(g2AA2 _ 1) 



(1.14) 



where Ai = t ~ ti, A2 



^(e2A(Ai+A2) _ 1) 

t, and Ai + A2 = ^2 - ti- 



2 Multivariate UOU-DifFusion Pricing Model 
2.1 Coupling UOU Processes 

Our goal is to construct a multi-asset price process (Sf)t>o with St = [S}, . . . , S*"), w here each individual 
asset price process {Sf)t>o, k = 1,2, . . . ,rt, is a univariate UOU diffusion obeying (1.4) with common 



drift parameter r and diffusion function a = . Suppose that each of the n univariate processes is 
described by its own set of positive parameters — {\k: J'fc; C fc, pfcj , fc = 1, 2, . . . , n. We denote by p| 



and a the univariate risk- neutral transition PDF of the form (1.10) and the diffusion coefficient given 



by (1.9), respectively, which both correspond to the fc-th asset price process. In the following, for each 
k — 1, 2, . . . , n, Up^, F'^, and X'^ will respectively denote the generating function, the mapping function, 

the inverse mapping function of the fcth diffusion model. The transition PDFs p\ and p^-^'"^^ correspond 

to the underlying diffusion {X^ )t>Q and the transformed diffusion (^X['""^^^ respectively AU the 

above functions indexed by k are obtained by using the parameters from S^k in the respective equations 
provided in Section [l] 

Recall that the processes {S^)t>o^ k — l,2,...,n, are defined by the transformation F*^' given by 



■'t 

(|1.8l), i.e. = F''{X^'"-''^) where {X^^'°'''')t>o are defined by (|l.5|) with X = Xk, v = Vk, and p = Pk 



For an infinitesimal time increment St, consider SS^ = S^_^_gf. — and SX^'^'"''^ = X^'^gf^ — xj.'""''\ It 

follows that 5S^ ~ {f''y{X^'""'''>)-SX^'""''\ From the map, (f'^Yix) = a''{f''{x))/vk, we readily see that 
the respective correlations between the ratios of the infinitesimal increment and the diffusion coefficient 
remain the same after the change of variables: 



6S^ SSi \ (SXI'"'-'''^ SX^"'-^^^ 



k,t = l,...,n. (2.1) 



This relation shows how to couple UOU diffusions. First of all, we can directly couple the processes 
^j^^(pfc,fe)^ ^ which, in turn, introduces correlations among the asset price processes {S^). In fact, the same 
method is used to couple geometric Brownian motions that are just exponentially transformed Brownian 
motions with drift. Hence, by using dependent Brownian motions, one introduces the correlations 
between the log-returns. In the case of more general diffusions, one may introduce correlations between 
the volatility-scaled returns. Alternatively, one may couple the underlying A"-processes. For example, it 
is possible to use a multivariate Ornstein-Uhlenbeck process as an underlying vector process. Another 
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approach (which is not discussed in this paper) is to extend the diffusion canonical transformation 
method directly to the case of multivariate diffusion processes. 

Another idea, an d the one that we implement in this paper, is the coupling of X'^''-' -bridges. As is 



shown in Subsection 1.3 the bridge PDF of a UOU X'^''' -process is the same as the density function 
of the Ornstein-Uhlenbeck bridge, which is a multivariate normal PDF. The respective multivariate 
distribution function obtained with the use of a Gaussian copula is nothing more than a multivariate 
normal CDF. Therefore, the coupling of X^'^^-bridges considerably simplifies the form of the joint path 
distribution function as well as the path simulation algorithm. 

2.2 Copula Function 

A copula C{ui,U2, ■ ■ ■ ,Un) is a multivariate cumulative distribution function (CDF) defined on the n- 
dimensional unit cube [0, 1]" such that every marginal distribution is uniform on the unit interval [0, 1] 
(for a more detailed definition we refer to [T5]). Suppose that <i>^, <i>^, . . . , $" are univariate distribu- 
tion functions, e.g., $*''(x) = f''{y) dy, where /'^ is a respective univariate PDF. It follows that 
C($^(a;i), . . . , $"(a;„)) is a multivariate CDF with marginals ^''{x) = C(l, .. .,l,Uk = ^''{x), 1, . . . , 1), 
k = 1,2, . . . ,n. The well-known Sklar's theorem states that any n-dimensional joint distribution function 
$ with continuous marginals . . . , $" has a unique copula representation: 

^ixi,X2, ...,Xn)^ C{<P\xi), $2(a:2), . . . , $"(a;„)). (2.2) 

In other words, for continuous multivariate distributions, the marginal distributions and the multivari- 
ate dependence structure can be separated. The multivariate density function f is then obtained by 



differentiating Eq. (2.2 1 



OUi . . . OUn 

As a corollary of Sklar's theorem we have a class of copula functions constructed from continuous 
multivariate probability distributions as follows: 

C{u^,U2,...,Ur.)^^{{'^^r\u^),{'i'^r\u2),...,{^n-\u^)) , (2.4) 

for any (mi, it2, . . . , u„) S [0, 1]", where ($'^)^^ is the inverse of . 

The Gaussian (or normal) copula is one of the most important in financial applications. This copula 
is constructed from the multivariate normal distribution: 

Cf'^''{u^,U2,. . . ,U^) = Mr{M~\ui),M~\u2), . . . ,M~\un)) , (2.5) 

where Mr is the standard n-variate normal CDF with mean vector zero, unit variances, and correlation 
matrix R. Here, A/""^ stands for the inverse of a standard univariate Gaussian CDF. 



Consider the problem of sampling from a multivariate distribution given by (2.2). The modelling 
algorithm consists of two steps. First, we simulate a uniformly distributed vector {Ui,U2, ■ ■ ■ ,Un) G 
[0, 1]" from the copulaC. After that, we sample a vector {Xi,X2, ■ ■ ■ , Xn) from the marginal distributions 
by evaluating the inverse CDFs: X^ — {^'')~^{Uk), k — 1,2, . . . ,n. As is seen, the crucial part of this 
algorithm is the efficient computation of the inverse of a distribution function. 
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2.3 Multivariate Path Copula 

The objective of this subsection is the derivation of a multivariate path distribution function in closed 
form. Such a function allows us to obtain an exact path simulation method and also to construct the joint 
path density function, which is used for calibrating the model with the maximum likelihood estimation 
method. 

We employ the copula method to construct the joint distribution function of the multivariate process 
X|''^ = {X^P''^\...,xj^''"''^^). The multi-asset price process St = {Sl,...S^) is then obtained by 

applying the respective mapping function to each univariate X'^^^-diffusion: S!^ = F'^' ^X^'''"'''''^ , k — 
l,2,...,n. 

Suppose that the process (X|''') conditional on Xq'''' is to be sampled at a set of times T = {tj}jLi G 

[0,T], T > 0, so that = to < ti < t2 < ■ ■ ■ < = T. Let T = represent some arrangement 

of time points in T. The ordering of the time points is determined by the simulation method used. 
For the (forward) sequential method we assume that 0<ii<---<tAr=T, i.e. Vj > ij — tj. 
For the backward-in-time bridge method we have that < <jv < t^^i < ■■■ < t2 < ti = i.e. 
Vj >^tj— tM+i-j- In other words, we first obtain the value of the process at the terminal time T and 
then simulate a bridge path conditionally on the previously sampled value and the initial value. Notice 
that for the bridge sampling method, the ordering of intermediate time points t2, ■ ■ ■ ,t]^ may vary, e.g. 
one can use a forward bridge sampling with < t2 < < ■ ■ ■ < tj^ < ti — T or a full bridge sampling 
method, where we successively halve the time interval or its largest segment to sample the process at 
the midpoint of the time segment conditionally on the previously sampled values at the endpoints of the 
current time (sub-)interval. The full bridge sampling may be applied along with the quasi-Monte Carlo 
method. 

Let fj{x) denote the PDF of X^^'''^''^^ conditional on the cr-algebra J^j^i generated by the previously 
sampled j path points x^'""''\ X^'""''\ . . . , X^'""'^-', where 1 < j < and 1 < k < n. For the sequential 



path sampling method, with ti = ti, we have that fj{x) — pj^*" — A^v"^ ' , x '^ 

For the backward bridge method we have that fi{x) = p'-^'"^'^ (T;Xjf'"'^\xj. Since all univariate 



processes are Markovian and are sampled backward in time, for each j ~ 2,3, . . . , N, fj{x) is a bridge 
PDF 6^ (ij-, x) of the Ornstcin-Uhlenbeck bridge conditional on x'lf'"'''' and x'^'''"'^\ where ij — <Ar+i-j 
and tj^i ~ As is shown in Subsection 1.3, t he PDF is a normal density function with 

kj ' 

A2 = tj^i — tj. 



mean fi^j and variance 6L-, with values given by (1.14| where A = A^, k — k^, Ai = tj — to ~ tj, and 



In A-space, the joint CDF ^j and respective joint PDF f,- of the n-dimensional point X^'''', j 



1,2, . . . , N , are then constructed by employing equations ( 2.2 ) and ( 2.3 1, respectively, where the marginal 

^ fk 
-oo •> j 



distribution functions are ^^{x) — fjix') dx' , k = 1,2, . . . ,n, j — 1,2, . . . , N , x e M., t > 0. Notice 



that, for the backward-in-time bridge method, the CDF is a normal CDF Af (^ ^ ^^.''^ ^ , for each 



From the Markov property of process (X^ ), the multivariate path distribution function <ll> and the 
respective joint path density function IF of ( X^''"' J conditional on Xq'''* are given by the respective 

\ / 7 = 1 N 
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products 



JV JV 

^(xi, . . . ,XAr) = J]^ *j(xj) and [F(xi, . . . , xat) = fj-(xj) 



(2.6) 



where Xj = (xj, . . . , a;"), j — 1,2 . . . ,N. In 5-space, the joint path PDF P and the joint PDFs of 
S( . , j = 1, 2 . . . , iV, are respectively given by 



JV 



P(Si 



,Sjvj 



n 



(x'>')) 



(2.7) 



where = (sj, . . . , s") 



To sample from the distribution function in (2.6) we need a fast and accurate algorithm for inverting 



univariate CDFs. Usually, for non-standard distributions such as the PDF given by formula (1.10 1, 
such algorithms are quite computationally intensive. Therefore, the sequential approach becomes rather 
unfeasible for the path simulation of the multi-asset price UOU process (St). 

As was mentioned above, we apply the Gaussi an c opula method to construct the underlying vector 



process xj''' with transition PDFs of the form (1.61 instead of direct coupling of the S-space asset 



price processes. This is then followed by applying a nonlinear transformation to obtain the asset prices 

To minimize the number of numerical inversions of CDFs, we employ a bridge simulation approach, 
which is followed by the coupling of the bridge distributions. The application of the Gaussian copula 
to bridge PDFs from the UOU diffusion family leads to a multivariate normal distribution. Indeed, if 

each CDF (say, the CDF of the fcth X'^) -bridge) is a normal CDF TV ('^) , then the multivariate 



CDF $ given by (2.2) with the Gaussian copula in (2.5) is just a multivariate normal distribution 
function Mb. (^^^^jT^' ' ' ' ' ^"tT^ " ) with mean vector (/ii, . . . ,/i„)^ and covariance matrix DRD, where 
D = diag(6i, . . . ,6„). 



2.4 Path Sampling with a Bridge Normal Copula 

Consider the following backward-in-time bridge sampling algorithm for the exact path-simulation of the 
multi-asset price process on a discrete partition T of the time-interval [0,T]. First, we generate the 
discrete-time random process (X|^')j=i_2,...,Ar conditionally on Xq''-' = (X^(S'g), . . . , X"(S'J)) , and then 

obtain sample values of (84^)^=1,2 JV by changing variables. Let us denote Xj = Xf'''"''^ and Sj = S^. 

/c = l,...,n, j = 0,l,...,iV. ' 

Step 1. Apply the inverse mapping functions to obtain the initial values: 

X^=X\S^), fc = l,2,...,n. 



Step 2. Sample the terminal-time value Xjy ~ {Xlf, . . . , A^Jv) ^oni the copula (2.2) by employing nu- 
merical inversion of the CDFs 



V^^'-^\t-xI,x)Ax, fc = l,2,...,n. 

-00 
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(i) Draw a normal vector (Zi, . . . , Zn) from the n-variate normal distribution function A/jj 
with mean vector zero and correlation matrix R. 

(ii) Obtain uniform variates Uk = Af{Zk), fc = 1, . . . , n. 

(iii) For each k = 1,. . . ,n obtain X'^ = (<f>'')^^(C/fe). 

Step 3. Sample Xj = (A|, . . . , A") for each j = A — 1, . . . , 1 by applying the bridge normal copula 
method as follows. 

(i) Draw a normal vector (Zi, . . . , Z„) from A/}j. 

(ii) For each k — 1, . . . ,n set X'^ = ^kj + bkjZk, where fikj and bf,j are given by (1.141 with 



respective parameters A = Xk, k — k^, Ai — tj — to, and A2 = tj+i — tj. 

the resulting ' 
at each time point: 



Step 4. Map the resulting values of the multivariate discrete-time process fxl^' j into the asset prices 



X"^ ^ S'; = f'^iXf), J = 1, . . . ,iV, fc = 1,2, . . . , 



n. 



As is seen from the algorithm, Step 1 involves the numerical inversion of n distribution functions. Since 
all parameters and initial spot prices are fixed, we can pre-compute and store the values of the inverse 
CDFs on a fine grid in (0, 1) and then apply a spline interpolation to sample X^r (see [13]). 

Notice that the exact bridge simulation method presented above is faster than any approximation 
method such as the Euler scheme or the Milstein scheme (e.g., see [H]). First of all, our method has 
no limitations on time increments which can be very large. Therefore, if a path needs to be sampled 
only at a few time moments, no intermediate times are required to guarantee the convergence of sample 
paths. Second, the Euler approximation method (or any similar one) being applied to an SDE defined by 



( 1.5 ) requires frequent computations of special functions that appear in the drift and diffusion coefficient 
functions. The bridge sampling method has only two computationally expensive steps, namely the 
sampling of terminal asset prices and the computation of mapping functions. As is mentioned above, the 
respective probability distributions can be tabulated to speed up the sampling. The mapping functions 
may be tabulated as well. 

The approach presented here can be applied to other families of hypergeometric diffusions (see [3 [8]). 
As is shown in [S], the probability distributions of the squared Bessel process and CIR process, which are 
the underlying diffusions for the so-called Bessel and confluent families of F-diffusions, respectively, are 
reduced to randomized gamma distributions. Such distributions are mixture probability distributions 
that can be obtained by allowing the shape parameter of the gamma distribution to be random. To 
couple F-diffusions, one can just couple the gamma distributions using a copula function. The use of a 
multivariate gamma distribution is more preferable, but all known examples of such a distribution admit 
only positive correlations. 

3 Calibration of the Multivariate UOU Model 

3.1 Univariate Case 



It is very important from the practical point of view to develop a reliable and reasonably quick calibration 
scheme for the UOU diffusion family. Our objective is to obtain a calibration scheme that provides two 
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levels of calibration: first, an initial full calibration of all parameters of the model and, second, a much 
faster recalibration that can be used as soon as new data have arrived. The second calibration scheme 
may be used throughout the day or even for longer periods, while the full calibration only need to be 
executed if markets move considerably. 

To estimate a best-fitted parameter set ^ = {A, i^, c, p} of the UOU model based on (observed) market 
option price data, the least squares method is employed. Suppose that an option with strike Ki and 
maturity Ti has an observed price Oi, while the model produces a price of Ci = C{Ki, T^; ^) for the same 
option, where i = 1, 2, . . . , M. The goal of the calibration process is to minimize the least squares error 
for the M options considered: 

M 

FiO ^yw^ \C{K,,Tf,0 - O.l' ^ min, (3.1) 



=1 



where Wi is a weight that reflects the relative importance of reproducing the ith option price precisely. 

The suitable choice of the weight factors Wi, i = 1,2, ... , M, is crucial for good calibration results. 
The confidence in individual data points is determined by the liquidity of the option. The weights can 
be evaluated from the bid-ask spreads: Wi = \Of^^ — 0^"^|~^. Alternatively, as it was suggested by [S], 
one may use the Black-Scholes (BS) "Vegas" evaluated at the implied volatilities of the market option 
prices to compute the weights: Wi = (^dC^^ {af^) / dcr^ , where dC^^/dcr denotes the derivative of the 
BS option pricing formula with respect to the volatility a, and af^ = a^^{Oi,Ki,Ti) is the BS implied 
volatility for the observed market price Oi. 

In general, the calibration of a pricing model is an inverse problem, whose solution depends discon- 
tinuously on the data. To achieve uniqueness and stability of the solution, a penalty function is added 
to the least squares term: 

M 

P^iO =^2"^^ \C{K„Tf,0 - O.I' +ai7(P,Po) ^ min, (3.2) 
1=1 ^ 

where the penalty function H is chosen such that the problem becomes well-posed. 

As is examined in the relative entropy method may be applied for solving ill-posed calibration 
problems. The relative entropy of a probability measure P on sample space fl with respect to some 
primal measure Pq is defined as follows: 



" dp ■ 




In 


-1 


dPo. 


Jn 



dP 

m — dP. (3.3) 



The regularization parameter a in (3.2) is used to adjust the trade-off between the accuracy of 
calibration and the numerical stability of results with respect to input option data. The right choice of 
a is based on the Morozov discrepancy principle [TU], which is described by the following algorithm: 

1. Compute parameters of of the primal measure Pq by solving the nonlinear least squares problem 



(3.1 ) in low precision. 



2. Fix S G (1, 1.5) and numerically solve equation -F'a(^o) = SFi£.o) for the regularization parameter 



a, where -F'q(Co) is defined in (3.2 1. 
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Figure 2: Market call price surface for IBM, July 7th, 2009 (left plot). Comparison of historical option 
prices and option prices calculated using the UOU model with the optimal parameter set (right plot). 



3.2 Numerical results for the univariate case 

The data set used consists of 79 European call option prices with maturities ranging from less than 
one month up to 1.56 years. These market prices were obtained from Yahoo for IBM having the spot 
share value of 101.34 on July 7th, 2009. For the sake of simplicity, the risk-free interest is assumed to 
be constant and equal to r = 0.25%, and the dividend rate is set to zero. The calibration routine was 
developed using Matlab with the Optimization Toolbox, running on an Intel Core 2 CPU 2.14GHz with 
2 GB of main memory. 

To obtain the set of parameters for the primal probability measure, the UOU model is calibrated to 
the historical data from May 7th to July 7th, 2009. Using historical asset prices, St-, j = 0, 1, . . . , TV, 
~ to <<i < ••• < tpf, and the transition densities, we obtain the following (single-asset) log-likelihood 
function for this set of observations: 

N 



Here, for simplicity, we assume the sequential simulation method. In case of a general sampling method, 
the log-likelihood function is given by 



(3.5) 



where fj is defined analogously to in Subsection 



2.3 



as if there was only one asset. 
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Table 1: Initial values and bounds for the parameters of the UOU model. 



Parameter 


P 


V 


c 




Lower bound 


0.001 


0.005 


45 


0.5 


Upper bound 


0.5 


2 


250 


10 


Initial value 


0.04 


0.34 


102.59 


1 



In practice, the implementation of the calibration procedure is started with some initial values of 
parameters. The upper and lower bounds for the parameters should also be provided. Based on the 
empirical analysis, such bounds are obtained and are provided in Table [T] 

The first step of the calibration procedure takes approximately 200 seconds to fit the model to 
63 historical asset prices. The optimal values, that maximize the log-likelihood function (3.4), are 
p — 0.0357, V — 0.0531, c = 118.2404, k, = 0.5951. This set of parameters defines the primal probability 
measure Pq. The estimation of the regularization parameter a is based on the algorithm described 
above. The calculated value of a is 0.266. 

The final step of the calibration process is the minimization of the nonlinear least squares function 
regularized by the relative entropy as is given in (3.2 1. The computation algorithm utilizes the Matlab 
function Isqnonlin with the exit tolerance set to 10"'^. This function employs the Levenberg-Marquardt 
least-squares algorithm for estimating optimal parameters. The starting values and the limits for the 
parameters remain the same as given in Table [T] The computational algorithm takes approximately 400 
seconds to fit the model to 79 option prices. The best-fitted parameters of the model are p = 0.0203, 
V = 0.0013, c = 102.1384, n — 0.6579. The objective function Fa attains its minimum value of 1.58. 

Notice that the discrepancy between the computed option prices and observed option prices may 
originate from different sources. First, the market data may contain errors or misleading information. 
For example, the values of illiquid options might be mispriced, or simple input errors may occur. Second, 
the calibration procedure estimates model parameters of an arbitrage-free model, while the market prices 
are not necessarily arbitrage-free. Hence, there is an inherent mismatch between the model prices and 
the market data. Notice that the use of time-dependent parameters may decrease the level and number 
of errors and make the calibration procedure maturity-wise. Another possible solution to improve the 
accuracy is to employ the calibration separately for out-of-the-money, at-the-money and in-the-money 
options. 



3.3 Multivariate Case 

Let us consider the multi-asset price processes (St)t>o modeled as described in Section [2j i.e. n UOU 
diffusions are coupled via the Gaussian copula function. 

The calibration procedure can be split into two stages: (1) estimation of the parameters of the 
marginal (single-asset price) processes; (2) estimation of the correlation matrix R of the Gaussian copula. 
Such a calibration algorithm admits multiple variations. First, one may use the maximum likelihood 
estimation (MLE) to fit the marginal models to historical asset prices. Second, one may use the least- 
square method to fit the marginal models to historical derivative prices (say European options). For 
both approaches, the correlation matrix is then estimated by the MLE using historical asset prices. 
Alternatively, one may use only observed asset prices to estimate all parameter of the multivariate 
model simultaneously without splitting the calibration process. Notice also that the multivariate path 
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distribution depends on the simulation method used. By using the sequential sampling or some version 
of the bridge sampling, one may obtain different models and, hence, obtain slightly different estimates 
of the model parameters. 

Let I S j = (Sf^^ , . . . , 5^ ) I be the n x (N +1) matrix containing + 1 independent historical 

prices for each of the n financial assets observed on a set of time points T = {to, ^i, . . . , tAr}. Let 
(^j^) = {£,1, ■ ■ ■ ,£,n, R) denote the set of parameters to be estimated. The historical observations in 
A^'^^-space are obtained by applying the inverse map: X^. = X''{Sf.;£k)- Suppose the joint path PDFs 



of the n-dimensional processes (xj'''') and (St) is constructed with the Gaussian copula as given by (|2.6[) 



and (2.7), respectively. The (n-asset) log-likelihood function is 

N 



L„(^,i?) =lnP(Si,...,SAr) = ^lnf;(S,) 

N 

= 5]ln0«(AA-i(<f>](Ai ;ei)), . . . M'\<f]{^;U}) (3.6) 



n N 



where denotes the joint PDF of the n-variate normal distribution with mean vector zero, unit 



variances, and correlation matrix R; Li is the single-asset log-likelihood function given by (3.4) or (3.5 1, 
and L'^" denotes the log-likelihood function for the copula function. Recall that the expression in (3.6) 
is independent of the simulation method used. For the sequential and bridge methods, we provide below 
specific expressions of the log-likelihood function. 



As is suggested by the structure of the log-likelihood function in (3.6), the calibration process can be 
split into two steps. First, the sets £k = {A/c, Vki Cfc, Pk} , fc = 1, 2, . . . , n of parameters of the marginal 
distributions are estimated by employing the maximum likelihood estimation: 



N 



argmaxLi(^fe) = argmax> In ^ J j i^i'^ £k) , fc = 1, . . . ,n. (3.7) 



As is seen, the parameters of the marginal distributions are estimated based on historical data. An 
alternative approach to computing the parameters is to fit asset price distributions to observed option 
prices as described in Subsection |3.1| 

The last step is the estimation of the correlation matrix R for the given optimal univariate model 

parameters ^ = {^i, . . . , ^fc = ^Xk, i^k, Ck, Pfe|, estimated during the previous step. 

First, we consider the sequential calibration method with the following log- likelihood function 



N 



'm)^Y.^nM^-\^]iXl;i,)),...,M-\<i>]{X-;U))- (3-8) 



Sequential Calibration. For the sequential path generation method, the algorithm is as follows. 
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(i) Map all the observations into X^''^ -space using the respective inverse maps: 

X, EE (Xj, ...,X^)^ ((Xi(§i.;6), • . . ,X"(^,'^.;en)), J = 0, ... ,7V. 

(ii) Compute vectors = (uj, . . . , u") e [0, 1]", = 1, . . . , A^, by evaluating the integrals: 



p^^^'=)(t,-t,_i;Aj=_i,x;a-)dx, fc = l,...,7i 



(iii) Maximize the log-likelihood function with respect to R: 

N 

LT'{R\l)^Y.^n(l>R{U-\u]),...,U-\u^))^m^^. 

The estimation of the log-likelihood function for the sequential calibration involves numerous estima- 
tions of the CDF for the UOU model. Since there is no simple-form solution for the CDF, the numerical 
integration of the probability density function should be performed regularly. By applying the bridge 
approach to the construction of the multivariate path distribution function, the number of integrals to 
be computed numerically on step (ii) reduces from n x iV to n. This is due to the fact that for the bridge 
approach, the CDF j = 2, . . . , iV, is Gaussian. Hence, for the bridge path generation method, the 
log-likelihood function can be simplified as follows: 

N 



J=2 



bkj bkj j (3.9) 

+ ln0«(AA-i(a>}(Ai 5 6)), . . . ,AA-i(<f?(X?;C„))), 
where mean /i^j and variance hf.- computed by formulae in (1.141; and ti — T. 

Bridge Calibration. The following algorithm can be applied for the backward-in-time bridge path 
generation method. 

(i) Map the observations into X^^^-space Sj — > Xj = F^^(Sj;^),j — 1,...,N, as is described in 
part (i) of the sequential algorithm. 

(ii) Compute un = (u^, . . . , u^), the values of normal CDFs corresponding to the terminal point of a 
path: 

u%= f " p^i''-''\tN-tN-i;X^N-i,x;^k)dx, k^l,...,n. (3.10) 



(iii) For each k = 1, . . . ,n and j = 1, . . . , N — 1 calculate Hkj and bkj by using ( 1.14 1 with respective 
parameters A = Aj., k — k^., Ai = tj — Iq, and A2 = ij+i — tj. Then, set 



^k ^ - ^^kJ 
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(iv) Maximize the log-likelihood function with respect to R: 



LT'm) - 5] ln0fl (a:], . . . , xj) + \ncbH {^-'{u],), . . . ,N-\u"n)) ^ max. 



Table 2: Optimal parameters estimated for IBM, Microsoft, Pepsi and Wallmart 





IBM 


Microsoft 


Pepsi 


Wallmart 


p 


0.0496 


0.2173 


0.0865 


0.0493 


V 


0.0887 


0.0365 


0.1149 


0.0886 


c 


103.9904 


21.1638 


31.671 


52.3842 


k 


0.9670 


0.874 


0.910 


0.9874 



3.4 Numerical results for the multivariate case 

For this numerical experiment the daily observations of four American companies, namely, IBM, Mi- 
crosoft, Pepsi, and Wallmart, have been collected from Yahoo!"'"'^. The examined period is April 7th, 
2009, to July 7th, 2009, and it consists of 63 time points. In the first stage of the calibration, the optimal 



sets of parameters of the marginal distributions are estimated by solving Eq. (3.7), and they are provided 
in Table H 

Two approaches are then used for the evaluation of the optimal correlation matrix R. In the first 
approach, the correlation matrix is obtained by the pairwise calculation of the correlation coefficients. 
There are (2) correlation coefficients for 4-asset price processes to be estimated. In other words, instead 
of so lving a 4-dimensional least-square problem with the log-likelihood function L'^°" given by (3.8) or 
(3.9), we independently solve (2) 2-dimensional problems with a 2-by-2 correlation matrix of the form 



and then compose another 4-by-4 correlation matrix using the estimations of 6. 

However, the resulting matrix may violate the positive-definite property. To overcome this problem, 
a method suggested by |14| of finding the closest correlation matrix by the spectral decomposition is 
applied. The resulting matrices R are shown in Table [s] The computation time is 32.6 seconds for the 
bridge simulation and 47.9 seconds for the sequential simulation. 



Table 3: Correlation matrices obtained by using the bridge path simulation (left table) and the sequential 
path simulation (right table). The pairwise computation of the correlation coefficients is employed. 



/ 1 
0.297 
0.151 
V 0.337 



0.297 
1 

0.089 
-0.045 



0.151 
0.089 
1 

-0.080 




/ 1 
0.278 
0.243 

V 0.336 



0.278 
1 

0.184 
-0.050 



0.243 
0.184 
1 

-0.051 



0.336 
-0.050 
-0.051 
1 



In the second method, the correlation matrix as a whole is estimated. The computation of an optimal 
correlation matrix is performed in Matlab using the function fmincon, which allows us to find a minimum 
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of a multivariate function with non-linear constraints. By adding nonlinear constraints, the algorithm 
works in the class of semi-positive matrices, which is absolutely necessary for the correct formulation of 
the correlation matrix. However, the candidate matrix, which minimizes the objective function in ( |3.6[ ), 
may not have ones on the principal diagonal. To obtain a correct correlation matrix that is closest to 
the given one, the spectral decomposition method is applied again. The results are shown in Table |4] 
and Table m 



Table 4: The candidate semi-positive matrix that minimizes the objective function L'j^" in (3.8) (left 
table) obtained by using the sequential simulation method and the closest correlation matrix obtained 
by using the spectral decomposition method (right table). 

/ 1.000 0.277 0.169 0.335 \ / 1 0.278 0.243 0.336 

0.277 0.993 0.127 -0.049 | spect. [ 0.278 1 0.183 -0.049 

0.169 0.127 0.480 -0.024 dc^p. 0.243 0.183 1 -0.035 

V 0.335 -0.049 -0.024 0.993 / \ 0.336 -0.049 -0.035 1 



Table 5: The candidate semi-positive matrix that minimizes the objective function L'^" in (3.9) (left 
table) obtained by using the bridge simulation method and the closest correlation matrix obtained by 
using the spectral decomposition method (right table). 



1.000 
0.290 
0.145 
0.335 



0.290 
0.973 
0.084 
-0.043 



0.145 
0.084 
0.966 
-0.076 



0.335 
-0.043 
-0.076 

0.993 



spect . 
decomp. 




0.293 
1 

0.086 
-0.043 



0.148 
0.086 
1 

-0.078 



0.336 
-0.043 
-0.078 
1 



4 Pricing Path-Dependent Options 
4.1 Path-Dependent Multi-asset Options 

Suppose the price process (St) is modelled at a discrete set of time points (written in an increasing 
order) T = {0 = ^0:^1:^2, ■ • ■ i^Af = T}. Hence, we construct a multivariate discrete-time n-dimensional 
price path (Si)i=i,2,...,Ar, where = (S'j, S^, . . . , S^). 

Let us consider two discrete-time monitored path-dependent securities, namely a Bermudan option 
and an Asian option. For an Asian-style option the payoff function Ajv(Si, . . . , S^v) is assumed to be 
a function of averages A'^ = A (S'f , . . . , 5^) of the asset prices, where, for example in the case of the 
arithmetic time-averaging, A'^ — X^i^i ^i- The Asian basket option with the payoff 

An = (mk^A''-K^ , (4.1) 

considered in Section [4. 2[ is an example of an arithmetic average option. 

The arbitrage-free value V — V{So,T) of a discrete-time monitored path-dependent option (without 
early exercise opportunities) can be represented as a mathematical expectation under a risk-neutral 
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probability measure Q: 

V = e-'^^EQ[A^(Si, . . . , Sm)] = e-'^^EQ[A^.(F(Xl), . . ■ ,F(X^))] . (4.2) 
To find the option value one may use the Monte Carlo method. 



The fair value of a Bermudan option cannot be represented in the form (4.2 ), but is a solution of some 



dynamic programming problem (see (4.3 1 below). Let A(t, S) denote the payoff function for exercise at 
time t in state S = (S"^, . . . , S"*) . One may consider the two following examples of the time-homogeneous 
payoff function A: 



• a max call option with the payoff A(S) = (niax(S'-'^, . . . , S"") — K 



+ 



a geometric average call with the payoff A(S) = (11 ^ ) " — K 

k=l 



Let V{t,S) denote the value of the option at time t given St = S, assuming the option has not 
previously been exercised. We are interested in the present value V{tQ,So). This value is determined 
through the backward-in-time recursion: 

V{tN,S) = A(i^,S) 

ViU,S) = max{A(t„S),E«3[e-(*'+i-*')y(t,+i,S,+i)|S, = S]}, ^^'^^ 

where i = N — I, . . . ,0. The risk- neutral expectation value function in the latter expression is called the 
continuation value function in state S at time 

In our computational tests, the univariate models correspond to the model parameter values used to 
plot the local volatility function with the thick line in Figure [T] So, for the UOU-family we simply set 
Pk = 0.02, Kk — 1, Vk — 0.5 and Ck = 100, for all k. The interest rate is r = 0.05. 

4.2 Pricing Asian Basket Options 

In this example we use the Monte Carlo method for pricing the arithmetic Asian basket call option with 



provides the results 



payoff function (4.1) under the multivariate UOU model. The Monte Carlo simulations were run on the 
SHARCNET network (http : //www. sharcnet . ca). Mostly we use Gulper — a 42-CPU cluster. The code 
was written in FORTRAN-90 using the MPI library. 

First, we present numerical results for a bivariate {n — 2) UOU model. The correlation matrix used 

in the normal (Gaussian) copula has the usual form of i? = ^ ^ 1 ) ' "^^^^"^ 

of our numerical tests for the two-asset model. We used the following values of parameters: number 
of observations N — 100, terminal time T — 1.0, number of scenarios M = 10 000 000, and spot value 
^0 = 100. 

The next example confirms that the computational complexity of the bridge copula method scales 
linearly with the number of assets. For all tests we use the crude Monte Carlo method on a single CPU. 
Clearly, the computational cost can be significantly reduced by using variance reduction methods and 
low-discrepancy point sets along with the parallel computing (see [7]). Another method of speeding up 
the computer code is to tabulate all distribution, generating and mapping functions that depend on the 
parabolic cylinder functions and other special functions. 

Table [7] provides the results of our numerical tests for increasing numbers of assets. We use the same 
values of parameters as those of the previous test but reduce the number of scenarios to A/ = 1 000 000. 
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Table 6: Pricing the Asian basket call option under the bivariate UOU model for different values of the 
correlation coefficient using the Monte Carlo method. A sample standard error is given after ± sign. 



K 


e = -0.75 


e = Q 


e = 0.75 


90 


31.538 ±0.008 


27.942 ±0.008 


22.444 ±0.009 


100 


22.786 ±0.007 


20.409 ±0.008 


16.168 ±0.007 


110 


15.614±0.007 


14.348 ±0.007 


11.321 ±0.006 



Table 7: Pricing the Asian basket call option under the n-variate UOU model using the Monte Carlo 
method. A sample standard error is given after ± sign. 



K 


71 = 2 


71 = 3 


n — b 


71 = 10 


90 
100 
110 


24.292 ±0.026 
17.627± 0.024 
12.403 ±0.021 


34.420 ±0.027 
25.974 ±0.026 
18.786 ±0.024 


45.504 ±0.028 
36.195 ±0.027 
27.507 ±0.026 


59.785 ±0.028 
50.274 ±0.028 
40.800 ±0.028 


Time 


3 745 sec 


5 340 sec 


8 389 sec 


15 307 sec 



A 10-by-lO correlation matrix R given in Figure [3] is generated using the random Gram matrix method 
from 12 . For values of ri < 10 we let R be the n x n-submatrix in the upper-left corner of the matrix 
given in Figure [3] Clearly, such a submatrix is again a correlation matrix of lower dimensionality. 



1.000 
0.550 

-0.123 
0.148 

-0.239 
0.238 
0.239 
0.096 
0.289 

-0.325 



0.550 

1.000 
-0.223 

0.188 
-0.109 

0.396 
-0.132 
-0.118 

0.237 
-0.390 



-0.123 

-0.223 
1.000 
0.198 
0.135 
0.242 
0.013 

-0.731 
0.244 

-0.134 



0.148 
0.188 
0.198 
1.000 
-0.519 
0.314 
-0.191 
-0.332 
0.469 
0.194 



-0.239 
-0.109 

0.135 
-0.519 

1.000 
-0.279 

0.292 
-0.044 
-0.077 
-0.443 



0.238 
0.396 
0.242 
0.314 

-0.279 
1.000 
0.023 
0.050 
0.049 

-0.062 



0.239 
-0.132 
0.013 
-0.191 
0.292 
0.023 
1.000 
0.423 
-0.306 
-0.292 



0.096 
-0.118 
-0.731 
-0.332 
-0.044 
0.050 
0.423 
1.000 
-0.417 
0.152 



0.289 
0.237 
0.244 
0.469 

-0.077 
0.049 

-0.306 

-0.417 
1.000 

-0.534 



-0.325 
-0.390 
-0.134 

0.194 
-0.443 
-0.062 
-0.292 

0.152 
-0.534 

1.000 



Figure 3: The 10-by-lO randomly generated correlation matrix R. 



The idea of jl2| is to generate independent pseudo-random vectors Ui, . . . , u„ distributed uniformly 
on the ?i-dimensional unit sphere and then to use the Gram matrix R = U^U, where U = (ui| • • • |u„) 
has Ufc as kth column and is the transpose of U. To create in R" we use a vector of independent 
standard normals zj. ^ Af{0,I) normalized by the Euclidian norm: U/j — zi;/\\z,f^\\. 



4.3 Pricing Bermudan Options 

Regression-based methods are broadly applied to pricing multivariate American options (see HT' and 
references therein). The main idea consists in the use of regression to estimate continuation values 
from simulated paths. Each continuation value is approximated by a linear combination of some basis 
functions. The coefficients of such a representation are estimated by a regression method (typically 
least-squares) . 
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Table 8: Comparison of price estimates for Bermudan and European put options on the maximum of 
two assets modeled by the bridge copula multivariate UOU model with different values of the correlation 
coefficient 6. A sample standard error is given after ± sign. 

9 Bermudan European 
-0.75 3.336 ±0.002 1.831 ±0.002 
7.106 ±0.003 5.680 ±0.003 
0.75 11.580 ±0.004 10.798 ±0.005 



Consider an expression for the continuation value of the form: 

L 

jeQ s,+i) I S, = S] = ^ = Pli^iS), (4.4) 

9=1 

for some basis functions ip = {ipi, . . . , '4'q ■ 1^" ~^ ^ and constants /3iq, q = 1,2,...,L, i = 

0, 1, . . . , A'^— 1. An approximate value Vq for the American (i.e. Bermudan) option price can be calculated 
by the following regression algorithm. 

Regression-Based Pricing Algorithm. 

(i) Simulate M independent asset price paths (Sy, . . . , Snj), j — 1, . . . , M. 

(ii) At terminal nodes, set VNj — A(tjv, ^Nj), j — ^, ■ ■ ■ , M. 

(iii) Apply backward-in-time induction for i = N — 1, . . . ,1: 

• Given estimated values Vi^i,j, j = 1, . . . , Af, the least-sc^uared estimate of /3i = (/3ii, . . . , PiR)^ 

is given by Pi = B^'^^B^v,i , where B^^i is the ixL matrix with (r, g)-entry ^ i^r{^ij)i'q{^ii) 

j=i 

and B^v,i is the L-vector with qih entry — — ^ %pq{Sij)Vi+i^j. 

i=i 

• Set % = max {^{U, S,^), Y.%q^q{S^j)^ , j = 1, . . . , M. 

(iv) Set = (Vn + . . . + Vui)/M. 

To illustrate the regression algorithm above, we consider a Bermudan put option on the maximum 
of two underlying assets 5*1 and ^2 modeled by the UOU bridge copula with interest rate r = 0.05, zero 
dividend yield, strike price K — 100, time to maturity T — 1.0, initial spot 5*0 — 100, and A^ = 10 
exercise times distributed evenly in the time interval [0,r]. Here we use the same choice of the model 
parameters as in the previous numerical example. Clearly, this example can easily be extended to the 
case with three or more assets. 

In the regression algorithm the basis functions are the power functions 1, 5*1, S2, Sf, S"!, 6*1 6*2, and 
payoff function A(S'i, 52) = {K — max(S'i, S'2))''^ . We apply the regression method with M = 100000 
independent paths. To calculate standard error we replicate the random estimator 100 times. Table [8] 
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shows numerical results for differing values of the correlation coefficient 9. For comparison's sake we also 
calculate the values of the European put option on the maximum of underlyings by the Monte Carlo 
method with 10^ sample paths. 

Conclusion 

In this paper we have constructed a new multi-asset pricing model, whose marginal (single-asset price) 
processes are nonlinear (local volatility smile) regular diffusions from a recently developed UOU family 
of probability conserving martingale models [51 H] . The multivariate model is based on a bridge copula 
method, where a normal distribution function couples the underlying univariate Ornstein-Uhlenbeck 
bridges and consequently forms a multivariate asset price process with built-in correlations. Such an 
approach preserves the solvability of the model, hence the multivariate path density is available in 
closed form and can be used for the calibration of the model to market prices. Extra flexibility of the 
multivariate model is provided by different variations of the bridge sampling method. We are also able to 
sample multivariate paths from their exact distribution. Moreover, the proposed exact bridge simulation 
algorithm runs faster than any approximation scheme (e.g., the Euler method). 

To illustrate the financial applications of our model, we succeeded in calibrating the model to single- 
asset equity option market prices as well as in calibrating the multi-asset price correlation matrix to 
historical asset prices. We also succeeded in pricing multi-asset Asian and Bermudan options by using a 
path integral Monte Carlo (MC) approach and an MC regression method, respectively. The preliminary 
calculations presented in this paper pave the way to further applications of derivatives pricing under 
such multi-asset state dependent diffusion models. The numerical efficiency of the path integral MC 
approach used in this paper can be further improved by employing quasi-MC methods. These methods 
have already been successfully used for another nonlinear model, the so-called Bessel-K model f5ll8], and 
are also applicable to the UOU local volatility smile model. 
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